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Transition to chaos in magnetized, weakly coupled plasmas 



(N 

o 



o 

a 

^— > 

a 

C 

o 
o 



Andrea Carati* 
Department of Mathematics, 
Universita degli Studi di Milano 
Via Saldini 50, 20133 Milano, Italy 

Francesco Benfenati'l' 
Corso di Laurea in Fisica, Universita degli Studi di Milano 
Via Celoria 12, 20133 Milano, Italy 

Alberto Maiocchi* and Luigi Galgani§ 
Department of Mathematics, 
Universita degli Studi di Milano 
Via Saldini 50, 20133 Milano, Italy 

Matteo Zuin 11 

Consorzio RFX, Associazione EURATOM-ENEA sulla Fusione, Padova, Italy 
(Dated: September 20, 2012) 

We report the results of numerical simulations for a model of a one component plasma (a system 
of N point electrons with mutual Coulomb interactions) in a uniform stationary magnetic field. We 
take N up to 512, with periodic boundary conditions, and macroscopic parameters corresponding 
to the weak coupling regime, with a coupling parameter T = 1/64. We find that a transition 
from order to chaos takes place when the density is increased or the field decreased so that the 
ratio LUp/uJc between plasma and cyclotron frequencies becomes of order 1 (or equivalently the ratio 
vi/Xd between Larmor radius and Debye length becomes of order 1). The result is in agreement 
with the theoretical prediction obtained in [1], on the basis of an old estimate of Iglesias, Lebowitz 
and MacGowan[2] for the intensity of the electric field acting on one electron and due to all the other 
ones. A comparison can be made with the threshold obtained from kinetic theory arguments, which 
corresponds to the condition v ee /u) c — 1, where v ee is the electron collision frequency. The latter 
threshold has a completely different dependence on the physical parameters and, for F = 1/64, gives 
a critical value of lj p about 80 times larger. 

PACS numbers: 05.45.-a, 52.55.-s 
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The microscopic foundations of plasma physics are usu- 
ally formulated in terms of kinetic theory, in which a key 
role is played by the concept of "collision frequency" or by 
the related one of "mean free path" . It is usually stated 
that for the so-called weakly coupled plasmas, such as 
gaseous-discharge plasmas, fusion plasmas, or a plasma 
in the solar corona, the Coulomb coupling is so small that 
"their thermodynamic properties are analogous to those 
of an ideal gas" (see [3], page 10), i.e., the electrons be- 
have essentially as if they were free. 

A different approach was taken in paper [1], in which 
the microscopic Newton equations themselves were tack- 
led directly, making use of the tools of ergodic theory and 
of a quite recent extension of Hamiltonian perturbation 
theory to the thermodynamic limit (see [4-6]). This al- 
lows one to obtain theoretical results for the microscopic 
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model itself, with no need of passing through the ap- 
proximation of the Boltzmann equation. In particular, 
for an infinite plasma immersed in a uniform stationary 
magnetic field B the electron motions were estimated to 
be ordered if the ratio lo p /lo c between electron plasma 
and cyclotron frequencies is below unity, chaotic in the 
opposite case. This means that the Coulomb interac- 
tions among the electrons are strong enough to produce 



a chaoticity threshold when lo p 



(the definitions of 



these and of some other familiar quantities will be re- 
called in a moment). Here "ordered" has to be under- 
stood in the sense of ergodic theory, i.e., that there exists 
at least one dynamical variable the time-autocorrelation 
of which does not decay to zero, or decays in an ex- 
tremely slow way. For completely chaotic motions, in- 
stead, the time-correlations of smooth dynamical vari- 
ables arc known to quickly decay to zero. 

Notice that the estimate for the chaoticity threshold 
determined in [1] can be eventually expressed in an ex- 
tremely intuitive way, namely, as the condition that the 
typical value of the perturbing force on any electron (the 
sum of the Coulomb forces due to all the other ones) just 
equals the typical value of the Lorentz force. So, denot- 
ing by Ej the modulus of the electric field acting on the 
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j-th electron and by Vj the modulus of the electron's 
velocity, the condition for the chaoticity threshold takes 
the simple form Ej ~ Bvj/c (in Gauss units). In turn, 
the electric field acting on one electron and due to all 
the other ones, looked at as a random variable, obviously 
has vanishing mean, so that its typical value is estimated 
by its standard deviation. The value of the latter, at 
density n e and temperature T, was estimated long ago 
by Iglesias, Lebowitz and MacGowan [2] to be given by 
y/An n e k B T, where ks is the Boltzmann constant. This 
leads for the threshold to the condition uj p /uj c ~ 1, which 
in particular is independent of the coupling parameter T. 

We used here the definitions u> c = eB/mc, u> p = 

\J e 2 n e /m 1 T = e 2 /aksT, where m and e are the electron 

mass and charge, c the speed of light, and a = n e 1//3 the 
mean interparticle distance. A relevant related quantity 
is the Debye length X D d = ^/ksT ] /n e e 2 . 

The theoretical estimate for the chaoticity threshold at 
uj p ~ u c found in [1] was quite unexpected because kinetic 
arguments apparently suggest that, in the weakly cou- 
pled regime r< 1, a transition might occur at v ee ~ lo c , 
where v ee is the electron collision frequency (see for exam- 
ple [7]). On the other hand one has (see for example [3], 
page 35) v ee ~ T 3 / 2 | logT| w p , and this gives a threshold 
at u p ~ co c / (r 3 / 2 | logr|), i.e., at u p > uj c (for T < 1). 

In this letter we report the results of numerical simu- 
lations at T — 1/64. A transition to chaotic motion is 
seen to occur for ui p /uj c in the interval between 0.25 and 
2, in agreement with the prediction given in [1]. 

Let us recall that a one component plasma model is 
just a system of N point electrons with mutual Coulomb 
interactions. We denote their position vectors by xj, j = 
1,...,N, and take them to lie in a box of side L (so 
that the electron density is given by n e = N/L 3 ). They 
are subject to the Lorentz force (e/c) B A Xj due to a 
constant homogeneous magnetic field B (which we take 
directed along the z axes), and to their mutual Coulomb 
forces. The electric force on the j-th electron, which 
depends on the positions xi, . . . , x^v of all electrons, may 
be simply denoted by e 2 E(x :) ), where E is the electric 
field acting on that electron and due unit charges located 
in the positions of the other electrons. As we are using 
periodic boundary conditions (so that we are actually 
dealing with a system of infinitely many electrons), the 
latter field can be computed [8] by the Ewald summation 
of the field due to an infinite cubic lattice of charges of the 
form Xi+Ln. Here, n is a vector with integer coordinates, 

i.e., n = (n x e x + n y e y + n x e z ), with n x , n y and n z e Z, 
while e x , • • • are the unit vectors along thr axes. 

Rescaling time by the electron cyclotron frequency lu c 
and position vectors by the mean interparticle distance a, 

i.e., introducing r = ui c t and yj = xj/a, the equations 
of motion take the form 



(the dots denoting now derivatives with respect to r), and 
so contain only one (dimensionless) parameter, namely, 
lu p /lu Ci while the rescaled density is obviously equal to 
1. In the simulations, the explicit form of the Ewald 
rcsummcd field E acting on the j-th particle is given by 

E (yj) ^mir^ls [erfc(ar;, n ) + exp(-a 2 rf n ) 

n I l n '"l L 

47r ^-^ ^-^ k k 2 . 

+ zs 2^ 2^ Tk]2 ex P(~i^ sm ( k ' r ') • 



def 



def 



y; =e z A Yj + E(y,) 



Here rj =' yj - y u while r Lxl =' y 3 - yi + Ln/a: the 
function erfc(x) is the usual error function, and a is the 

Ewald convergence parameter which we chose as a d = 
7r 1 / 2 jV 1 / 6 _L _1 . In the first sum the term corresponding 
to the self-force on the j-th particle should be excluded. 

These are the equations of motion that were actu- 
ally integrated numerically, using a symplectic splitting 
method. The conservation of energy in every run was 
better than a part over 10 3 . The integration time was 
chosen proportional to uj c in order that all different cases 
be integrated for the same physical time. In any case, 
the time was always some hundreds cyclotron periods. 

The initial data were chosen in the following way: the 
electron positions y^- were taken uniformly distributed in 
the box of side N 1 / 3 , while the velocities were extracted 
from a Maxwellian with a given temperature T. This in- 
troduces in the model (in addition to uj p /uj c ) the further 
parameter T, or equivalently the dimensionless Coulomb 
coupling parameter T, to be used in the Maxwell distri- 
bution for the velocities. 

For what concerns the number N of electrons in the 
box, our computational power allows us to go up to 
N = 512. This induces a lower bound on T, namely, 
r > N~ 2 / 3 . Indeed, in order to correctly simulate the 
Coulomb cumulative force acting on an electron in a 
plasma, the side of the box has to be at least equal to 
the Debye length, which, in our rescaled units, takes the 
value \ D = V- 1 / 2 . We took T = N- 2 ' 3 . Computations 
were performed both for N — 128 and N — 512, which 
correspond to T = 128~ 2 / 3 ~ 0.04 and T = 512~ 2 / 3 = 
1/64 ~ 0.016 respectively. 

We now come to the main issue, i.e., whether the 
motions are ordered or chaotic. Obviously what plays 
the role of the unperturbed system with completely or- 
dered motions is the limit case with u> p /u) c — 0, for 
which the Coulomb interaction disappears and one has 
pure Larmor gyrations. The problem then is to deter- 
mine whether a threshold for chaotic motions takes place 
as the parameter uj p /uj c is increased and T is varied. 
To this end we considered the magnetization of a box, 

M = (e/2mc) J^Xj A xj, looking at its autocorrelation 
function (normalized by NksT) 

dcf < M(0)M(t) > 

CM{t) = Nk^f ' 
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FIG. 1: Top: Autocorrelation Cm(£) of magnetization versus time for oj p /uj c = 0.25 (left) and u) p /u) c = 2 (right). The time 
scale is the same in both figures, and one has u> c = 1 at the right, cj c = 8 at the left. Notice the fast decay to zero at the right. 
Bottom: Discrete Fourier transform (absolute value) of C.m(£) versus lo/uj c for lo p /lo c = 0.25 (left), and oj p /u c = 2 (right). 
Peaks (and thus also magnetization) have disappeared at the right. Here T = 1/64. 



and at its Fourier transform Cm (to) . The latter is a phys- 
ically very relevant quantity because, according to lin- 
ear response theory (see [9, 10], or Appendix B of [11]), 
iwCx(w) gives the susceptibility x( w ) at frequency u. In 
the formula for the time-autocorrelation Cm(£), the av- 
erage < • > is meant as a phase-average with respect to 
Gibbs measure; in our computations, however, we esti- 
mated it by the time-average along an orbit (with initial 
data extracted as previously explained), as often done in 



numerical works. We did not investigate the relations 
between the two averages. Moreover, the Fourier trans- 
form Cx( w ) was estimated by the amplitude of the dis- 
crete Fourier transform of Cm{^), which will be simply 
called the spectrum. So we report figures of the time- 
autocorrelation Cm versus t, and of the corresponding 
spectrum versus angular frequency uj/uj c . 

Having fixed T — 1/64, by increasing uj p /uj c we found 
that a threshold occurs for uj p /uj c between 0.25 and 2. 
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This is exhibited in Fig. 1, where the results are reported 
for such two values of u) p /u! c , 0.25 on the left and 2 on 
the right. The autocorrelations are reported in the upper 
part of the figure, and the spectra in the lower part. 

For cjp/uj c = 0.25 the autocorrelation is seen to dis- 
play regular oscillations with a decreasing amplitude: we 
were unable to follow this relaxation process up to the 
end. The oscillations are apparently peaked about the 
cyclotron frequency and its low harmonics (as should be, 
due to the nonlincarities in the equations of motions). 
This is clearly exhibited by the spectrum, with its large 
peak at lu/lo c = 1, and the smaller ones about the low 
harmonics lu/uj c = 2,3,.... Of special relevance is the 
peak at uj = 0, which corresponds to the existence of a 
nonvanishing static susceptibility, i.e., to the existence of 
diamagnetism. There also appears a continuous compo- 
nent, which accounts for the extremely slow drift towards 
equilibrium. This case clearly corresponds to prevalently 
ordered motions with a corresponding nonvanishing dia- 
magnetism, and should be interpreted as an indication 
that the perturbation due to the Coulomb interactions 
is not yet sufficiently large to produce prevalent chaotic 
motions. The passage to chaos, however, already oc- 
curred at ujp/uje — 2. Indeed in this case the autocor- 
relation is seen to go to zero in an extremely short lapse 
of time (even shorter than one cyclotron period 2it/uj c ), 
so that the peaks disappear from the spectrum and one 
only remains with the continuous part. This means that 
for r = 1/64 the threshold in u! p /u> c lies between 0.25 
and 2. For L = 128~ 2 / 3 the corresponding figures at 
those same values of uo p /uo c are qualitatively similar to 
the above ones, and are not reported here. 

So, the numerical results obtained for T = 128~ 2 / 3 ~ 
0.04 and L = 1/64 ~ 0.016 are in rather good agree- 
ment with the theoretical prediction found in [1], namely: 
at ujp/uj c = 1 the interactions become strong enough as 
to make the motions chaotic. On the other hand this 
is apparently in contrast with kinetic theory arguments, 
according to which Coulomb interactions should be neg- 
ligible up to values of uj p larger by a factor 20 and 80 



respectively. The discrepancy would become enormous 
in physically relevant cases, as gaseous-discharge plas- 
mas, fusion plasmas, or a plasma in the solar corona, for 
which T takes the typical values 10~ 3 , 10~ 5 , 10~ 7 respec- 
tively. Indeed, in terms of densities, for fusion plasmas 
the coupling should be negligible up to densities about 
thirteen orders of magnitudes larger than according to 
the law ujp/ujc = 1- Notice by the way that, as shown 
in [1], the latter threshold appears to fit pretty well, at 
least as orders of magnitude are concerned, the empirical 
data for disruptions in fusion machines. 

In conclusion, the present numerical work confirms, for 
weakly coupled plasmas, the theoretical predictions for 
the chaoticity threshold given in [1]. The main point is 
however that this confirms the estimate given in [2] for 
the intensity of the electric field acting on one electron 
and due to all the other ones, which turns out to be much 
larger than usually assumed. Now, the fact that the col- 
lective Coulomb effects are relevant for weakly coupled 
plasmas (with r <C 1) is very well known (see for ex- 
ample [12], page 8). Indeed it is just for such plasmas 
that the number r~ 3 / 2 of effectively interacting electrons, 
(those contained in a Debye sphere) turns out to be very 
large. What is apparently lacking, perhaps, is a general 
acquaintance with how large such an effect may actually 
be, in fact so huge as to possibly explain the disruptions 
in fusion machines. 

Such an acquaintance might perhaps help elucidating 
also the situation met in the problem of anomalous trans- 
port, where it occurs that "... measured energy transport 
rates typically exceed those calculated for binary collisions 
..." [13], or even "... greatly exceed" them [14]. We hope 
to come back to this point in the future. For a study 
on anomalous diffusion in a strongly coupled one com- 
ponent plasma, through numerical computations of the 
same type as those performed here, see [15]. 

Acknowledgments. The present paper is dedicated 
to Francesco Guerra (La Sapienza University at Rome) 
on the occasion of his seventieth birthday. 



[1] A. Carati, M. Zuin, A. Maiocchi, M. Marino, E. Martinez 
L. Galgani, Chaos 22, 033124 (2012). 

[2] C.A. Iglesias, J.L. Lebowitz, D. MacGowan, Phys. Rev. 
A 28, 1667 (1983). 

[3] S. Ichimaru, Plasma Physics: an Introduction to Statisti- 
cal Physics of Charged Particles, Benjamin (Menla Park, 
1955). 

[4] A. Carati, J. Stat. Phys. 128, 1057 (2007). 
[5] A. Maiocchi, A. Carati, Commun. Math. Phys. 297, 427 
(2010). 

[6] A. Carati, A. Maiocchi, Commun. Math. Phys. 314 , 129 
(2012). 

[7] R.D. Hazeltine, J.D. Meiss, Plasma Confinement, 

Addison- Wesley (Redwood City, 1991). 
[8] P. Gibbon, G. Sutmann, in Quantum Simulation of Com- 



plex Many-Body Systems: from Theory to Algorithms, 
J. Grotendorst, D. Marx, A. Muramatsu eds., NIC Se- 
ries 10, 467-506 (2002). 
[9] R. Kubo, J. Phys. Soc. Japan 12, 570 (1957). 

[10] Yu.L. Klimontovich, Statistical Physics, Harwood Aca- 
demic (Chur, 1982). 

[11] F. Benfenati, A. Carati, L. Galgani, Chaos 21, 023134 
(2011). 

[12] T.J.M. Boyd, J.J. Sanderson, The Physics of Plasmas, 

Cambridge U.P (Cambridge, 2003). 
[13] E.J. Doyle et al., Nucl. Fusion 47, S18 (2007). 
[14] J.W. Connor, HR. Wilson, Plasma Phys. Control. Fusion 

36, 719 (1994). 

[15] T. Ott, M. Bonitz, Phys. Rev. Lett. 107, 135003 (2011). 



